Diffusive high-temperature transport in the one-dimensional Hubbard model 
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We consider charge and spin transport in the one-dimensional Hubbard model at infinite tem- 
perature, half-filling and zero magnetization. Implementing matrix-product-operator simulations of 
the non-equilibrium steady states of boundary-driven open Hubbard chains for up to 100 sites we 
find clear evidence of diffusive transport for any (non-zero and finite) value of the interaction U. 
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I. INTRODUCTION 

The one-dimensional (Id) Hubbard model is the sim- 
plest model of strongly interacting dynamics of spinful 
fermions on a lattice within a single-band approximation. 
Its Hamiltonian reads 



H 



L-l 



E (4, 



*=i «e{t,4-} 



L 

,Ci+i )8 + h.c.) + U^n^nji, 

1=1 



(1) 



where the operators Ci, 8 , q s , for spin s e {1\4-} an d site 
i G {1, . . . , L} are standard fermionic annihilation and 
creation operators and rii )S :— c\ s Ci_ s are the density op- 
erators. The model is solvable by a coordinate Bcthc 
ansatz 1 and possesses an infinite number of conservation 
laws^. While stationary properties of the Id Hubbard 
model are well understood, see the monograph^, much 
less is known about its dynamics, for instance about the 
transport behavior. In thermodynamically the most in- 
teresting regime, at half-filled band and zero magnetiza- 
tion J2i=i( n i,s) — L/2, studied in the present paper, the 
model is gapped for charge excitations and gapless for 
spin excitations, for any U ^ 0. At zero temperature it 
is therefore an example of a Mott (charge) insulator and 
ballistic (ideal) spin conductor. Transport can be qualita- 
tively characterized by a Drude weight - a linear response 
indicator of ballistic transport, defined as the weight of 
zero-frequency singular term S(u>) in the real part of con- 
ductivity cr(cu). Spin and charge Drude weights at zero 
temperature have been calculated in Ref^, with finite- 
size corrections given in Ref£, while a regular part of 
a(uj) is studied in Rcf. 6 . At nonzero temperature on the 
other hand no rigorous result is known and there is no 
consensus between numerical and Bethe ansatz based re- 
sults. Thermodynamic Bethe ansatz suggests 7 - that, even 
at half-filling, the charge Drude weight is finite, so the 
model was predicted to exhibit ideal charge transport; 
for a similar conclusion see also Quantum Monte Carlo 
calculation in RefA Analytical calculations at large U on 
the other hand support vanishing charge Drude weight 9 , 
for a study of low-energy excitations see also Ref>!£. Ex- 
act numerical simulations of small systems, again at half- 
filling and at high/infinite temperature, suggest^ ~ 1/L 
scaling of the charge Drude weight. For temperatures 
much smaller than the gap semiclassical arguments to- 



gether with field-theoretical scattering rate predicts dif- 
fusive transport^. Vanishing finite-temperature Drude 
weight in thermodynamic limit (TL) L — > oo offers the 
possibility of an insulating or diffusive behavior. 

However, at high or infinite temperatures, non- 
equilibrium transport properties of Id Hubbard model in 
either charge or spin sector are not known as there has 
been up to date no analytical or numerical method capa- 
ble of reliably treating this regime. In this paper we em- 
ploy non-equilibrium steady state simulations^ using an 
efficient matrix product ansatz^ for the time-dependent 
density matrix and provide a clear evidence of diffusive 
transport for both attractive U < and repulsive U > 
cases at infinite temperature. Namely, we show 1/L scal- 
ing of charge as well as of spin current and clear linear 
density profiles. 



II. BOUNDARY DRIVEN HUBBARD CHAIN 

Using the Jordan- Wigner transformation we can map 
the Id Hubbard model ([T]) to a spin-1/2 ladder system. 

>W _ 



Namely, writing Cj-j- = P\VxC! i where P\" 



p(°\> p( r ) _- 



where P} 



spin-up fermions, and Cjj, 
rf • • ■ rf for spin-down fermions, one can verify that 
fermionic operators Cj jS , c\ satisfy canonical anticommu- 
tation relations provided af and r" are two sets of Pauli 
matrices (and af := (of ± iaJ)/2, if := (r* ± irJ)/2). 
Writing the Hubbard Hamiltonian (JTJ) in spin-ladder form 
one obtains 

L-l 

h = - ■= Y. (°l<+i + <%°Ux + ^i + T M+i) + 



u 



1=1 
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The spin-1/2 ladder system consists of two XX chains in 
two legs and a Z — Z type interchain coupling along the 
rungs. For numerical simulations of the Hubbard model 
we shall use this ladder formulation ©. 

To induce a noncquilibrium situation two legs are cou- 
pled to independent reservoirs. Their action is decribed 
in an effective way via the Lindblad equation^ for the 



density matrix p of the ladder system, 
d 



At 



p = i\p,B\ + C*>(p), 



(3) 



where the dissipative term is expressed in terms of Lind- 
blad operators Lk, as 



C dis ( P )=J2([L k p,Ll} + [L k ,pL{] 



(4) 



We use eight Lindblad operators acting locally on the first 
and last sites of each leg, injecting or absorbing fermions 
(spinons) with certain probability: 



for the first, and 



L 5fi = ^T(lTn)Tf, Lr t8 = v/r(l ± p) t 



L • 



(6) 



for the second leg. Y is the strength of the coupling 
to the baths while p is a driving strength playing the 
role of a chemical potential bias. As demonstrated in 
previous studies of Id spin chains^ the precise form of 
Lindblad operators does not influence the bulk proper- 
ties. Because of dissipative terms the time-dependent 
solution p{t) of the Lindblad equation converges after a 
long time to a time-independent state called a nonequi- 
librium steady state (NESS), p^ — lim t ^ 00 p(t), which 
is unique^. Once the steady state is reached, expressed 
in terms of a matrix product operator of a given bond 
dimension (following the method described in detail in 
ReLi2 straightforwardly adapted for the spin ladder), ex- 
pectation values of arbitrary observables in the NESS can 
be efficiently evaluated. All expectation values consid- 
ered in this paper are taken with respect to the NESS, 
that is (A) = tr (paoA), which we will - when it is clear 
from the context, and to simplify notation - denote just 
by A. In each NESS calculation we have carefully checked 
that the convergence is reached, i.e., we evolve the Lind- 
blad equation ([3]) until a time- independent state is ob- 
tained, and that the results are stable with respect to 
increasing bond dimension^. For p, = 0, i.e., no driv- 
ing, one has an equilibrium setting, resulting in a trivial 
NESS poo oc 1. This means that for small driving p we 
are studying nonequilibrium behavior at an infinite tem- 
perature. Note that such an infinite temperature state is 
separable in the operator space. Since the efficiency of 
the numerical method crucially depends on the entangle- 
ment infinite-temperature nonequilibrium states are the 
easiest ones to calculate because the entanglement is ex- 
pected to be smaller than at finite temperatures. 

Expectation values of fermionic observables are ob- 
tained from the corresponding ones in the ladder formu- 
lation, for instance, particle densities are n^ — (of + 1 )/2 
and riii = ( T i + -0/2- Magnetization currents of the two 
spin species, defined through the continuity equations 
d(a!/2)/dt = j! ct) - j£\, d(7f /2)/dt = j\ T) j£\, are 



Ji — 2\ a i a i+l a i a i+l)iJi ~ 2 \~i T i+1 T i T i+l)- 

In fermionic picture the particle (charge) current is j\ = 

jf +Ji , while the spin current is j\ = (jj — j^')/2. 
Particle density is m — n^ + n^, while spin density is 
Si = (riif — riii)/2. Because of the same driving at both 
ladder legs the currents j( a > and J*- 7- ) are the same. There- 
fore, NESS is such that it has a nonzero charge current 
and zero spin current. We have also performed simu- 
lations with Lindblad operators on the r-chain driving 
transport in the opposite direction, that is with 



£ 5 ,6 = \/r(l±M)Tf, L 7 ,8 = vr(i T li) T 



(7) 



In such a case of spin driving the NESS has a nonzero 
spin current and zero charge current because p T > — 
—p 17 ' holds. Furthermore, we stress that spin and 
charge transport are interchanged under the particle-hole 
transformation for the down spin fermions only and si- 
multaneously chaging the sign of U. Namely, taking 

R ■= IliLiTf = Rf one finds RJi S) R^ = 3i C} /2 and 
RH(U)B) = H(—U), provided one takes a symmetric 
interaction term (n^ — h)(nn — h) in ([T]) or, equivalently, 
adds a chemical term —UN/2 to H with N = ^ m <s - 
Even though our master equation evolution ([3]) does not 
strictly conserve N, we have checked explicitly that the 
results based on Hamiltonians H and H—UN/2 are iden- 
tical. 



III. RESULTS 

A. Evidence of diffusion: density profiles and 
scaling of currents 

We set t = 1, T = 1 and p = 0.2, except in Fig. [7] 
where p = 1. Driving p = 0.2 corresponds to equilib- 
rium density in the reservoirs of «l,s = 0.4 at the left 
end and rip{. s = 0.6 at the right end. The average filling 

ratio is therefore n = 1/2, X^»=i n *t = Si=i n 4 = L/2. 
The value p = 0.2 is at the upper end of a linear response 
regime. For large drivings p > 0.6 one gets a negative 
differential conductance effect, where the current de- 
creases with increasing driving. The main goal of this 
paper is to classify spin and charge transport, whether it 
is ballistic, diffusive or anomalous. For NESS, different 
transport regimes are reflected in the scaling of the cur- 
rent on the system size. Fixing the driving strength p, 
in a ballistic conductor the current is independent of the 
system length, j ~ L°, for a diffusive conductor it scales 
as j ~ 1/L, whereas in the anomalous case the current 
is proportional to a fractional power of L. We therefore 
calculated NESSs for different sizes L. Typical density 
profile is shown in Fig. [TJ One can see that the densities 
of spin-up and spin-down fermions are linear in the bulk. 
Jumps in the density at the boundaries are due to over- 
simplified Lindblad operators that are not "matched" to 
the bulk dynamics, i.e., there are boundary resistances. 
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FIG. 1. (Color online) Density profile m,a along the chain for 
L = 100, U = 1. Apart from jumps at the boundary, density 
is linear which is typical for diffusive conductors. Solid black 
line, overlapping with the numerical points, is a best-fitting 
linear function. 
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FIG. 3. (Color online) (a) Density profiles n; ia at U — 2, and 
(b) the corresponding connected density-density correlation 
function (ni-|-nij.) c . Data are shown for L = 16 (blue squares), 
L = 32 (green circles) and L = 64 (red triangles), all for 
charge driving. Inset in (a): scaling of the jump between the 
reservoir and the 1st particle and between the 1st and 2nd 
particles, with size L. Black dashed lines indicate ~ 1/L and 
~ 1/L ' 7 . Note that at L = 64 boundary jumps still account 
for around 25% of the total density difference between the 
chain ends. 



FIG. 2. (Color online) Scaling of charge current y' divided 
by the extrapolated density drop LVn with the system size 
L for different interactions U. Thick full (red) line, overlap- 
ping with U — 1.0 data, is ~ 30.4/L, indicating a diffusive 
transport. Two dashed lines also suggest ~ 1/L scaling. 



Because these jumps are rather large we have fitted a 
linear function to the density profile in the bulk, thereby 
obtaining the density gradient Wn^ = V»4 = Vrij/2. 
In Fig. [5] we then plot the scaling of the charge current 
(which is in the NESS independent of the site) with the 
gradient of the charge density. At interaction strength 
U = 1 one can see a nice scaling j^ ~ 1/L. Together 
with a linear density profiles this is a clear indication of 
diffusive charge transport. As mentioned, for spin trans- 
port virtually the same behavior is obtained (data not 
shown). For U = 2 the scaling is not quite as good. It 
seems that for shorter chains j( c ' decreases with L slower 
than 1/L, however, for two largest sizes that we managed 
to calculate, a crossover to ~ 1/L scaling is clearly sug- 
gested. For smaller interaction U = 0.5 the convergence 
seems better, but contrary to the U = 2 case, the current 
approaches the asymptotic scaling (~ 1/L) from above, 



i.e. it decays a bit faster for short chains. It is not clear 
whether U = 1 corresponds to a crossover point between 
the two behaviors since details of convergence in the ther- 
modynamic limit might depend on a particular choice of 
boundary Lindblad operators. In Fig. HJa) we show den- 
sity profiles for U — 0.5 which, apart from considerable 
boundary jumps, again look linear. 

Looking at the density profiles at U = 0.5, 1,2 which 
are linear in the bulk already for rather small sizes L, it 
seems natural to conjecture that the transport is diffusive 
in TL L — > oo for all nonzero finite values of U. Different 
transient scaling of the current with L for shorter chains 
is likely due to rather strong boundary effects. That the 
boundary effects are notable can also be seen in the inset 
of Fig. [UJa), where we show the jump in the density be- 
tween the reservoir and the first site nif — riLf, as well as 
between the first two sites in the system n,2-\ — n±-f. While 
the boundary effects show a tendency to disappear in TL, 
at U = 2 and largest length L = 64 they are still non- 
negligible. One can try to optimize the coupling constant 
r in order to minimize the boundary effects, however we 
found that r f=a 1 is usually close to the optimal value 
which does not seem to depend on L. For U = 0.5 the 
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FIG. 4. (Color online) Density profiles (a) and connected 
correlations (b) for U = 0.5. Other parameters are the same 
as in Fig. [3] 



boundary effects are larger than for larger U cases stud- 
ied, U — \ and U = 2, see Fig. HJa). This is probably 
due to a smaller bulk resistivity which makes the effects 
of the contacts (contact resistivity) relatively larger. 



B. Density-density correlations 

In Figs. [3^b), QJb), and also EJb), we show 
the connected spin-up spin-down correlation function 
(n,tn,|)c = (n^n.i) - (n^}{nii}, that gives on-site cor- 
relations between two fermion species. If we extrapolate 
our finite-L data to TL, we find 

(n it nu) c ex ( M 2 - {2{n it ) - l)(2(n 4 ) - I)), (8) 

with a proportionality prefactor depending on U only, 
yielding a parabolic correlation profile for our linear den- 
sity profiles. Interestingly, in the middle of the chain the 
connected correlations become independent of the system 
size, while they are going to zero at the boundaries^. In 
Fig. [5] we in addition show also, for U = 1 case, the 
non-connected correlations (ro^n^) (top frame (a), red 
squares), or centered non-connected correlations ((n^ — 
l/2)(n ti - 1/2)) (bottom frame (b), red dotted line). We 
can see that the non-connected correlations have similar 
shapes as density profiles. 

Considering spin driving ([7]), we find numerically the 
same diffusion constant as for charge driving ((6]), while 
the connected correlations (rii-fnii) c change sign (see 
Fig. EJb), blue circles). Note that this implies that con- 
ductivities at infinite temperature do not depend on the 
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FIG. 5. (Color online) Density-density correlation functions 
for the case of charge driving (red squares) and for spin driving 
(blue circles). Frame (b) shows the connected correlations 
(symbols), frame (a) the non-connected ones. In (b) we also 
show ((riif — \){nii — |)) (dashed curves). All is for U — 1, 
L = 32. 



sign of interaction U . We have also checked explicitly, 
by comparing data for U = — 1 with 17 = 1, that den- 
sity profiles, currents and correlations are insensitive to 
the sign of U. One can see that in the presence of spin 
current without charge current, non-connected correla- 
tions, shown in Fig.[5ja) with blue circles, are practically 
independent of the site and are slightly smaller than 1/4. 
Note that in both cases, of charge and spin driving, 
the connected correlations scale as ~ /j 2 and are of purely 
nonequilibrium origin, i.e. they vanish in the equilibrium 
limit (p = 0). 



C. Large interaction U 

In the limit U — > oo the low energy excitations of the 
half-filled Hubbard model can be effectively described 
by the Id isotropic Heisenberg model. In our open sys- 
tem formulation this mapping cannot be strictly imple- 
mented, due to the presence of high-temperature baths 
which drive the system locally away from half-filling. It 
is therefore an interesting question whether the transport 
properties of the Hubbard model in the limit of large U 
are qualitatively the same as for the Heisenberg model. 

In Fig. [5] we plot density profiles in our open Hub- 
bard model for increasing values of U, keeping L fixed, 
and find increasingly arcsin-like shape, similar as in the 
isotropic Heisenberg model which displays an anoma- 
lous transport 2 ^ with the magnetization current scaling 



as ~ \j\fL. This perhaps explains slower decay of the 




FIG. 6. (Color online) As one increases U at fixed length 
L = 32 the rescaled density profiles (2rii, s — l)//i, shown for 
U — 1 and U — 5, become similar to ^ arcsin [ 2% £ — 1] (full 
red curve), found in the isotropic Heisenberg model 2 ^. 

current with L in the Hubbard model for small L's and 
larger U, seen for instance in Fig. [5] at U — 2. Note that 
the limits U — > oo and L — > oo do not commute. In order 
to recover the Heisenberg behavior in TL one has to first 

let U — > oo and only then L — > oo. 



D. Strong driving, fi = 1 

In previous subsections we have shown that the weakly 
driven Hubbard model displays diffusive behavior. Here 
we show that for strong driving, where the system is far 
away from equilibrium, the behavior of physical observ- 
ables can be dramatically different. For example, we 
briefly discuss the case of maximal driving fi = 1 and 
find that the current scales sub-diffusively as j c oc 1/L 2 . 
The corresponding density profile is shown in Fig. [7] We 
can see that the profile is in TL given by a simple cosine 
shape n^s — sin 2 (-7r(2i — 1)/4L), exactly the same as has 
been found analytically in the isotropic Heisenberg model 
at strong driving 23 . This suggests that similar exact so- 
lution for NESS at maximum driving fj, = 1 as for the 
Heisenberg spin chain is also achievable for the Hubbard 
model and points to wider applicability of the algebraic 
method proposed in Refs i 23 ' 24 . 



E. Non-half-filled case 

So far, all the results shown were for symmetric driv- 
ing, producing on average half-filled bands. However, if 
the filling of the two fermion species is not 1/2 we ex- 
pect the transport to be ballistic at high temperatures. 
This follows from an existence of a nontrivial constant of 
motion, which in the non-half filled case possesses non- 
vanishing overlap with the spin/charge currents^. In 
order to numerically verify the consistency of our non- 
equilibrium setup with this expectation, we choose a 
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FIG. 7. (Color online) At maximal driving, fi = 1, density 
profiles have a cosine shape (red squares, for L = 50, (7 = 1) 
while the current scales as ~ 1/L (data not shown), exactly 
as in the Heisenberg model at maximal driving 23 . In the inset 
we show convergence of An = m, 3 — sin (nx/2) with L. 
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FIG. 8. (Color online) Density profiles (a) and currents (b) 
for a non-half-filled Hubbard system. We show data for L = 
16 (squares) and L — 32 (circles), all for U — 1. Density 
profiles in the bulk are flat and currents do not depend on 
L, indicating ballistic transport. In both frames upper two 
sets of symbols are for one fermion species, lower two for the 
other. 



nonsymmetric driving with L1.2 = \/r(l =F MLt) a i an d 
£3,4 = \/r(l ± //Rf) <J L for the first chain, where /U14- = 
0.5 and ^ Rt = 0.9, while L 5fi = y/T(l =F Mu) rf and 
L 7:8 = y/T{l ± /xrj.) t l with /UL4. = 0.4 and ^ RJ , = 0.8 
for the 2nd chain. Density profiles can be seen in Fig. |U 
One can see that the gradient is very small (or zero); 
what is more, currents are almost independent of system 



size L. Namely, in Fig. [5Jd currents are within numerical 
errors the same for L = 16 and L = 32 (small inho- 
mogeneities visible in the Figure are due to truncation 
errors). If the transport were diffusive, as is the case 
for the half-filled system, the current for L = 32 should 
be half as large as for L = 16. We therefore confirm 
that the non-equilibrium transport is clearly ballistic for 
a non-half-fillcd Hubbard model. 



IV. CONCLUSION 

Summarizing our findings about the transport in half- 
filled zero-magnetization Id Hubbard model at an infinite 
temperature, we have shown that at finite interaction U 
both charge and spin transport are diffusive. This con- 
clusion is based on the scaling of the currents with the 
system size for up to 100 particles, as well as on per- 
fectly linear density profiles away from the boundaries. 
This complements ballistic spin transport at U = and 
arbitrary temperature and anomalous transport at high 
temperature and U — oo. We acknowledge support by 
the grants Pl-0044 and Jl-2208 of Slovenian Research 
Agency (ARRS). 



1 E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968); 
Physica A 321, 1 (2003). 

2 B. S. Shastry, Phys. Rev. Lett. 56, 1529 (1986). 

3 F. H. L. Essler, H. Frahm, F. Gohmann, A. Kliimper, 
and V. E. Korepin, The One- Dimensional Hubbard Model 
(Cambridge, 2005). 

4 B. S. Shastry and B. Sutherland, Phys. Rev. Lett. 65, 243 
(1990). 

5 C. A. Stafford, A. J. Millis, and B. S. Shastry, 
Phys. Rev. B 43, 660 (1991). 

6 E. Jeckelmann, F. Gebhard, and F. H. L. Essler, 
Phys. Rev. Lett. 85, 3910 (2000); R. M. Fye, M. J. Nartins, 
D. J. Scalapino, J. Wagner, and W. Hanke, Phys. Rev. B 
44, 6909 (1991). 

7 S. Fujimoto and N. Kawakami, J. Phys. A 31, 465 (1998). 

8 S. Kirchner, H. G. Evertz, and W. Hanke, Phys. Rev. B 
59, 1825 (1999). 

9 N. M. R. Peres, R. G. Dias, P. D. Sacramento, and 
J. M. P. Carmelo, Phys. Rev. B 61, 5169 (2000). 

10 S-J. Gu, N. M. Peres, and J. M. P. Carmelo, J. Phys.: 
Condens. Matter 19, 506203 (2007). 

11 P. Prelovsek, S. El Shawish, X. Zotos, and M. Long, Phys. 
Rev. B 70, 205129 (2004). 

12 S. Sachdev and K. Damle, Phys. Rev. Lett. 78, 943 (1997); 
K. Damle and S. Sachdev, Phys. Rev. Lett. 95, 187201 
(2005). 

13 T. Prosen and M. Znidaric, J. Stat. Mech. 2009, P02035 
(2009). 

14 G. Vidal, Phys. Rev. Lett. 91, 147902 (2003); 
F. Verstraete, J. J. Garcia-Ripoll, and J. I. Cirac, 
Phys. Rev. Lett. 93, 207204 (2004); A. J. Daley, C. Kol- 
lath, U. Schollwock, and G. Vidal, J. Stat. Mech. 2004, 
P04005 (2004). 



15 V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, 
J. Math. Phys. 17, 821 (1976); G. Lindblad, Comm. Math. 
Phys. 48, 119 (1976). 

16 NESS is in our case unique and is therefore reached irre- 
spective of the initial state p(0). Evans theorem 1 ^ states 
that NESS (fixed point of a Liouvillian dynamical semi- 
group p)l) is unique iff the set of operators {H, Li,L,2, . . .}, 
in our case {H, <Tj l , t^ l }, generates, under multiplica- 
tion and addition, the entire algebra of operators, in our 
case over the spin ladder on L sites. Indeed, this can be 
verified by operator identities giving a recursive genera- 
tion of af: o+ = (l/4)al[a+,[H,al]], and for j > 2, 
<r+ +1 = -<?+_! - (l/2)a*[a-,cr+Ha+]. Operators for the 
second leg r are generated through a similar recursion, 
replacing a" by rf. Finally, all possible products of {a- } 
and {t, } generate the entire operator algebra of the Hub- 
bard chain. 

17 D. E. Evans, Comm. Math. Phys. 54, 293 (1977). 

18 We use bond dimensions up to D — 80. Convergence times 
toa typically scale oc L, for example for L — 100 we needed 
to simulate ((3]) up to t^ « 500. 

19 G. Benenti, G. Casati, T. Prosen, D. Rossini, and M. 
Znidaric, Phys. Rev. B 80, 035110 (2009). 

20 X. Zotos, F. Naef, and P. Prelovsek, Phys. Rev. B 55, 
11029 (1997). 

21 Note that small noisy irregularities in the parabolic corre- 
lation profile in Fig. [4jb) for L = 64, being of order 10 -4 , 
are due to insufficient matrix product bond dimension D. 
Nevertheless, we have checked that the steady state cur- 
rent, which is here j c (a 0.05, has already converged with 
respect to increasing D. 

22 M. Znidaric, Phys. Rev. Lett. 106, 220601 (2011); 
M. Znidaric, J. Stat. Mech. 2011, P12008 (2011). 

23 T. Prosen, Phys. Rev. Lett. 107, 137201 (2011). 

24 T. Prosen, Phys. Rev. Lett. 106, 217206 (2011). 



